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The Martian Climate: Energy Balance Models with CQ2/H20 Atmospheres 
Progress Summary 

The last semi-annual status report (March 1, 198-1 - September 1, 1984) 
described the necessity to consider the meridional dependence of the 
Martian seasonal pressure wave. The basis is that atmospheric mass will 
be transported by the same eddy diffusive mechanisms as the transport of 
heat. This report presents a technique for the investigation, and develops 
coupled equations for mass and heat transport in a seasonal Mars model 
with condensation and sublimation of C02 at the polar caps Up to the 
present, we have used the heat transport model to calculate latent-heat 
C02-mass fluxes between ground and atmosphere with seasonal and 
latitude resolution. These fluxes will be taken into the mass transport 
model in the next phase of the program The contents of this report are: 

(a) Physical Considerations of Planetary Mass and Energy Balance 

(b) Effects of Phase Changes at the Surface on Mass and Heat Flux 

(c) Atmos )heric Transport and Governing Equations 

(d) Numerical Analysis 

Seasonal Mars Model: Formulation for a Volatile Planetary Surface with 
Atmospheric Transport of Mass and Heat Parameterized bv Fddv Diffusion. 

(a) Physical Considerations of Planetary Mass and Energy Balance 

In formulating an energy balance climate model lor Mars we characterize 
a vertical zonally-averaged column composed of both the carbon dioxide 
atm: sphere and a thermally-interactive surface regolith layer by a single 
temperature T(t,$), where <t> is the latitude and t is time The atmosphere 
has an approximately exponentially decreasing density distribution p(t,$,z) 
- Pa(t/t>,z)e _z/h with altitude z from the surface at z * 0 to z -> <x>, where 
p a is the surface density and h the scale height; while the thermally- 
interactive surface extends from z * 0 down to some thermal pentration 
depth z = d, where d = (2 ks/^) 1/2 , where Ks is the regolith (soil) thermal 
diffusivity and 2tr/n is the period of the annual heating wave (Hoffert and 
Storch, 1979). 

The atmospheric part exerts a surface pressure p(<f>,t) in general. The 
equation of state in the atmosphere is p = p a RT, where R is the gas constant 
of CO2 Let: 
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m (4>,t) 5 1 pdz a f^h * p/g (1) 

Jo 

be the zonally -averaged mass of the atmospheric column, where the 
expression on the far RHS is from hydrostatic equilibrium, where g is the 
graviational acceleration at Mars' surface. The corresponding atmospheric 
scale height is h * p/(p a g) * RT/g. If Ca and Cs are the constant-pressure 
specific heats per unit mass of the atmosphere and regolith layers, the heat 
capacities per unit area of the atmosphere and ’soil- per unit surface area 
area, respectively, mc a = p a c a h and Cs = psCsd, respectively. 

(b ) Effects of Phase Changes at the Surface on Mass and Heat Flux 

Let L 2 RT 2 dfcnp/dT (Clausius -Clapeyron equation) be the latent heat of 
sublimation /condensation — approximately constant for Mars conditions — 
so the saturation vapor pressure over dry ice at the seasonal caps is 

p(T) = (constant) x e _L/RT . 

Over ice-free surfaces, the pressure is unconstrained by the underlying 
surface temperature, but is affected by transport from the seasonal caps 
which are continually subliming and condensing in response to seasonal 
insolation cycles. We assume CC>2-ice covered surfaces are in vapor 
pressure equilibrium with the atmosphere such that the ground 
temperature is 

Tg = T c = - [ (R / L )Hn (p /constant) 1 ' 1 * 150 K (1) 

whenever the surface is frost-covered. In particular, we set Tg = T when T 
> T c and Tg = T c when T < T c We also assume bare ground albedos a& and 
longwave emissivities e when T >.T C and and different ice-covered values 
04 and 6i when T < T c If S($,t) is the solar flux per unit area of the planet, 
and parameterizing the longwave cooling to space by a grey body 
approximation in terms of the ground temperature, we write the net 
radiative heating as 

S(<t>.t)t 1 - <x b ] - £ b ori; if T > T c 

F*(T,<t>,t) • (2) 

S(4>.t)l 1 - <Xj] - ej<TT c 4; if T < T c 

where 0 the Stefan -Boltzmann constant. Since ground temperature is held 


*■» 

- J - 


at T c over ice, a sublimation mass flux m* is injected locally into the 
atmosphere over ice to balance radiative and meridional eddy diffusive 
heating (regolith desorption and adsorption of carbon dioxide by processes 
other than sublimation/freezing are neglected) Assuming the sublimation 
heat flux Lm* is supplied by the net radiative heating over ice F*(T) plus a 
F diff diffusive meridional heat flux (see equation 6) gives 

0; if T > T c 

m*(T,4>,t) = (3) 

IS(4>,t)l 1 - ttjl - E[OT c 4 ♦ Fatirl/L; if T = T c 

where m* = am/dt = - dm^/at, in which mj is the C02-ice present on the 
ground. Notice that m* can be either positive (sublimation) or negative 
(condensation), depending on whether solar heat or infrared cooling 
dominates F*(T) over ice. The term F<jjff will be zero at points over the ice 
cap, and wiU be positive only at the ice edge; also, it is small everywhere 
on Mars compared to the solar heat and infrared cooling. Here we couple 
the atmospheric temperature to the ice, so that T is never less than T c . 

(c ) Atmospheric Transport and Governin g Eq uations 

Meridional transport of both atmospheric mass and heat is by some wind 
field V ♦ V, where the “primed” component represents "eddy" fluctuations 
relative to the unprimed zonal-, altitude- and short period time-average 
velocity field. The eddies give rise to corresponding fluctuations in column 
mass and temperature which vary about the means, such that local values 
are m ♦ m‘ and T + T*. By definition <V‘>, <m‘> and <T‘> are zero, but the eddy 
mass and heat fluxes <V‘m'> and <VT> are not negligible, and may 
dominate the effects of the mean meridional flow for complex 
three-dimensional atmospheric circulations. Notice that both mass and heat 
fluxes are carried by the same velocity fields, and hence by the same eddy 
circulations, whose effects are subsumed under the < >-average 

Integrating over the column and < > -averaging the results then gives 
columnar conservation of mass and energy equations in the form: 

dm 

— + V*[Vm + <V'm‘>] = m*(4>,t) 
dt 



8T 

lmc a ♦ Csl — ♦ mc a V*l VT ♦ <VT'>] ♦ Lm*(4>,t) = F*($,t) 
at 

The divergence operator for a 2onally -averaged planetary atmosphere 
treated as a thin sphered shell is: 

6[cos4»( )$] 

V»()s 

acos<frd$ 

where a is the planetary radius. 

It is likely that after vertical, horizontal and time averaging, the residual 
mean meridional velocity in the atmosphere is small compared to the eddy 
velocities, V « (V' 2 ), as is the case with the earth's general circulation. This 
suggests a parameterization of the eddy heat fluxes in terms of an eddy 
diffusivity k as in North et al (I960), James and North (1962), etc. 
However, and this is the critical new feature of this model, we assume that 
since atmospheric mass is transported by the same eddies as heat, that it 
has the same meridional eddy diffusivity, i.e.. 


- <V'm’> - <VT> 

K = = 

8V/a04> dV/ad4> 


Substituting these approximations, and defining a new meridional variable 
x = sin<$>, (4) 


proportional to the planetary surface area measured from the equator, 
gives (after some algebra), the conservation equations, 

6m 1 T 0 2 m 0m 1 

— = — j ( 1 - x 2 ) - 2x — j + m*(T,x,t), (5) 

0 t r m [ 0x 2 Qx J 

0T l f 6 2 T 0Tl Lm*(T,x,t) F*(T,x,t) 

— = — j ( 1 - x 2 ) - 2x — j + , (6) 

0t tt [ 0 x 2 0x J C C 


where C = mc a + C s is the heat capacity per unit surface area of the 
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atmosphere plus thermally -interactive regolith layer, and 

r m 5 a 2 /n, 1:7 h (C/mc a )a 2 /x (7) 

are timescales for meridional transport of atmospheric mass and heat 
Notice rj scales by an additional factor proportional to the ratio of 
atmospheric heat capacity to total heat capacity, which makes poleward 
heat transport relatively less important when the surface’s heat capacity is 
large relative the the atmosphere’s. With m*(T,x,t) and F*(T,x,t) given by 
equations (2) and (3) based on the local T(x,t), the simultaneous solution of 
equations (6) and (7) by explicit finite -differences is straightforward. 

Calculating m* during condensation or sublimation of C02-ice is 
accomplished as follows: If T = T c , then first assume m* = 0, and calculate 
6T/8t from equation (6). If this assumption results in 6T/6t < 0, then 
condensation of C02-ice is occurring. If, instead, 0T/8t > 0, and, 
importantly, m^ > 0, then C02-ice is subliming. In either case, the latent 
heat transfer holds the temperature constant. This procedure only 
determines that m* is not zero, but does not calculate its value; this is 
accomplished by setting 0T/0t = 0, and using equation (6) to calculate m*. 
The only other possibility if T = T c is that m* = 0, and 0T/0t > 0 after 
assuming m* = 0 in equation (6); then, in fact, m* = 0 and equation (6) is 
used to calculate 8T/0t. 



- 6 - 


(cl) Numerical Analysis 

Relevant parameter values for present Mars conditions as given by Hoffert 
et ai. (1981): 

a - planetary radius * 3.39 x 10& m 

p 0 ■ planetary mean s’trface pressure * 700 Pa = 700 kg/s 2 -m 

g - surface graviation acceleration * 3.73 m/s 2 

mo = planetary mean atmospheric mass/unit area * 159 kg/m 2 

c a = atmospheric heat capacity /unit mass * 830 J/kg-K 

c a mo = atmospheric heat capacity /unit area * 1.56 x 10 5 J/m 2 -K 

The heat capacity per unit surface area of the ground estimated from the 
thermal inertia fit to temperature histories at the Mars lander sites is 

Cs = regolith heat capacity /unit area * 1.3 x 10& J/m 2_ K 

The system heat capacity /unit area is therefore: 


C = Cs ♦ c a mo * 1.5 x 10& J/m 2 -K 


And the atmosphere /surface heat capacity ratio is 

(caimo/C^aLrs or 0. 1 

This contrasts with the atmosphere /land heat capacity ratio of ^ 6 on earth 
where the atmosphere is much thicker; but since the high heat capacity of 
the oceanic mixed layer dominates the net thermal inertia in the terrestrial 
case, the overall atmosphere /system heat capacity ratio on earth is (Hoffert 
etal., I960): 


(c a m 0 /C) ear th ^ 0.04, 


which is actually less than present-day mars. 

Stone's (1974) atmospheric baroclinic instability model applied to mars 
conditions gives a meridional diffusion coefficient (Hoffert et al , 196 1) D = 
(mc a /a 2 k ^ 0.027 W/m 2_ K, and a corresponding meridional eddy 
diffusivity k or a 2 D/m 0 c a ~ 2 x 10& m 2 /s. James and North (1982) suggest 

the lower value D * 0.004 W/m 2 -K, corresponding to k ~ 3 x 10 5 m 2 /s. 
For our present model we assume k ~ lx 10& m 2 /s for both heat and 
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mass transport, which is intermediate to these values. 

The corresponding timescales in equations (5) and (6) are: 

1 mars year 

r m = a 2 /K ^ l.lx 10 7 s x - 0.2 mars years 

59 x 10 7 s 

1 mars year 

Ty s (C/mc a ) a 2 /n ^ 1.1 x 10 s s x ~ 2 mars years 

59 x 10 7 s 

Notice that despite the same merdional diffusivities, meridional mass 
transport is some 10 times taster than meridional heat transport because 

the high thermal inertial of the regolith makes atmospheric transport 

relatively less important. Nevertheless, unlike James and North (1962) we 
do not assume this transport is infinitely fast (T m -mx>), since it occurs over 
a significant fraction of the martian year We show below that it also occurs 
over a time which is long compared with the timestep with which we must 
integrate the conservation equations, wtech indicates that equation (5) as 
well as equation (6) must be solved to find the seasonal pressure wave 
p(x,t) = m(x,t)/g. 

Equations (5) and (6) can be written in terms of a generic variatle q(x,t) as 

dq 1 ’ 0 2 q 6q" 

— =— (1 -x2) 2x — ♦ q*(T,x,t), 

9t Tq [ 9x2 

equal to m(x,t) or T(x,t), where Tq is the relevant timescale for meridional 
mixing and q*(T,x,t) is the relevant mass or heat source per unit area. In a 
finite-difference approximation, we assume the domain from x = - 1 to x = 1 
is discretized into K intervals of width Ax = 2 /K, such that the Jth point is 
located at x(J) = -1 + 2J/K, where J = 0, 1,2,.. ,K , and J = 0, K are the endpoints 
at the south and north poles, respectively. Let I be an index denoting the 
time t = I At, where At is the timestep, from some inital time t = I = 0 when 
the array q(0,J) is known. The object is to find q(I,J) at all future times 
Assuming a time-explict diffrencing scheme and centered differences then 
gives q(I+ 1,J) at all points except the endpoints I = 0 and I = K in the form: 

At 

q(M,J) = q(I,J) ♦ •lA(J)q(I,J+l) ♦ B(J)q(IJ) ♦ C(J)q(IJ-l)l ♦ Atq*(IJ), 

TqAX 2 
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where 

A(J) = [1 -x(J) 2 - x(J)Ax], 

B(J) = 2|X(J)2 - 11, 

C(J) = [1 - x(]) 2 * x(J)Axl. 

The endpoints at J = 0 and J = K are evaluated by recognizing that second 
derivatives drop out since they are multiplied by [1 - x(J) 2 ], and by 
evaluating the first q derivative by one-sided differences: 

q(I+ 1,0) = q(I,0) ♦ [At/rq]((2[q( 1 ) - q(0)]/Ax}l + At q*(I,0) for x = -1 

q(M^) = q(I,K) - [At/r q ][{2[q(K) - q(K-l)]/Ax}] + At q*(I,K) for x = -1 


The stability criteria for explicit solution of such parabolic PDEs is the 
Courant condition. At < i/2TqAx 2 . For a spatial stepsize of Ax = 0.1 (K = 20 
points) we therefore require At < 0.0004 mars years for equation (5) and 
At < 0.004 mais years for equation (6). These timesieps are clearly much 
less than those required for mass (pressure) relaxation. Thus we need some 
250C computational steps/mars year to find the m(x,t) field but only 250 
steps/mars year for the T(x,t) field. However, the mass field relaxes to 
repeating cycles over a time of 0.2 mars yr, while the temperature relaxes 
in about 2 mars years. This suggests a strategy in which the temperature 
field is computed first — by running (say) 5 annual cycles to be 
conservative - and the source term array m*(I,J) from the last year is 
saved. The weak coupling of the m -equation to the T -equation through the 
dependence of -cj on the local m(x,t) may be neglected in view of 
uncertainties in k by replacing m by the longterm mean value m 0 . (In fact 
both r m and rj can probably be held constant in the calculation.) The mass 
field can then be computed with an order of magnitude smaller time 
resolution over (say) two annual cycles, saving only the last one 

Results can be plotted as contours of constant T and constant p = m/g on 
the t,x plane, or (L s ,x) plane where Ls is the aerocentric longitude, where 
intersections with the latitudes, or x*s, of the VL- 1 and VL-2 landers give 
the time-histories of T(t) and p(t) at the lander sites. To show the seasonal 
forcing, it would be nice to see contours of constant insolation S on the 
(L s ,x) plane as well. Finally, the pressure histories from the lander 
barometers at the two Northerm hemisphere locations should be 
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cross -plotted against observations and compared with the r m -♦ 0 
approximation of James and North (1962). 
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